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We present a new scheme for constructing initial data for the Einstein field equations using the 
conformal thin-sandwich formulation that does not assume conformal flatness or approximate Killing 
vectors. This includes a method for determining free data based on superposition, as well as a way 
to handle black hole singularities without excision. We numerically solve the constraint equations 
using a multigrid algorithm with mesh refinement. We demonstrate the efficacy of the method with 
initial data solutions for several applications: a quasi-circular binary black hole merger, a dynamical 
capture black hole-neutron star merger, and an ultrarelativistic collision. 



I. INTRODUCTION 

The purview of numerical relativity has extended to 
include not only relativity theory, but a wide range of 
other topics. Motivated by current and upcoming ef- 
forts to detect gravitational waves [1-5] , there has been 
extensive work on mergers of binary compact objects 
(BCOs) [6] including binary black holes (BH-BH) [7-9], 
binary neutron stars (NS-NS) [10] and black hole-neutron 
star (BH-NS) systems [11, 12]. In addition to binaries in 
quasi-circular orbits there has also been recent studies 
of eccentric binaries as may arise from dynamical cap- 
ture [13-15]. Other work of interest to astrophysics in- 
clude gravitational collapse of stars [16, 17], black hole 
accretion [18], and the nature of cosmological singulari- 
ties [19, 20]. Aside from astrophysical systems, numerical 
relativity has also emerged as a useful tool to explore vari- 
ous concepts in gravity and high energy physics [21], such 
as critical collapse [22] , ultrarelativistic collisions [23-27] , 
the gauge/gravity duality [28-32] , gravity and black holes 
in higher dimensions [33-35], and the (in)stability of Anti 
de-Sitter spacetime [36]. In all these applications, a nec- 
essary ingredient is a good method for constructing ini- 
tial data (ID). Here we present a new initial data solver, 
based on the conformal thin-sandwich (CTS) [37] for- 
mulation, which we have designed to be more generally 
applicable to a range of physical scenarios by avoiding 
symmetry or simplifying assumptions. 

There has been extensive research on the problem of 
constructing ID for general relativity and detailed re- 
views can be found in [38-41]. Early attempts at solving 
the initial data problem relied on certain assumptions to 
make the mathematical formulation of the problem more 
tractable, such as conformal flatness and maximal slic- 
ing. The widely used Bowen-York solution [42] is one 
such example. These assumptions are restrictive since, 
for example, the isolated Kerr black hole docs not admit 
conformally flat slices [43] and consequently the Bowen- 
York solution cannot be used to construct black holes 
with spin higher than S/M\ DM = 0.928 [44]. Other ex- 
amples include the use of quasi-equilibrium assumptions 
for constructing ID for binary systems (such as approxi- 
mate helical Killing vectors or the like, and approximate 



hydrostatic equilibrium for any matter in the system) ; see 
for example [45-60]. This serves as a good approxima- 
tion for astrophysically motivated quasi-circular inspiral 
sufficiently far from merger, though is not valid for eccen- 
tric mergers or the ultrarelativistic scattering problem. 
Many of these studies made further simplifying assump- 
tions, such as conformal flatness, which does not have an 
astrophysical motivation. Attempts to supply more real- 
istic conformal initial data include superposition of iso- 
lated black hole spacetimes [61-65], and in addition using 
post-Newtonian solutions [66, 67] and matched asymp- 
totic expansions to supply an initial outgoing radiation 
field [68-70] . Using superposed data allowed evolution of 
binary black holes in quasi-circular orbits with spins ex- 
ceeding the Bowen-York limit [71] . A further alternative 
approach, initially applied to binaries including neutron 
stars, involves solving the full Einstcin-Euler system of 
equations with a waveless and/or near-zone helical sym- 
metry approximation [59, 72-74]. 

Since our goal is to have a more general purpose numer- 
ical initial data solver that can be used for a range of ap- 
plications, as outlined in the first paragraph, we will use 
the CTS formalism with arbitrary conformal metric and 
other free data to be chosen as needed for the particular 
application. For our first version of the code, as presented 
here, we restrict to four-dimensional, asymptotically-flat 
spacetimes, with application to BCO interactions. For 
the free data we use superposed, boosted single CO 
spacetimes. At large separation this is a good approx- 
imation for the physical metric of dynamical capture bi- 
naries and the ultrarelativistic scattering problem, and 
the non-linear corrections from solving the CTS equa- 
tions are small. For quasi-circular binaries, again at large 
separation this is a good approximation. However, un- 
like the scattering problems, at practical (due to limited 
computational resources) initial separations to allow evo- 
lution through merger, the simple superposition we use at 
present will not give improved astrophysically relevant ID 
compared to current quasi-equilibrium approaches. Com- 
pared to existing studies using superposed data, what is 
novel about our work is the inclusion of matter to also 
model CO fluid "stars," and the consideration of ultra- 
relativistic initial boosts with Lorentz factors up to 10. 
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Another novel aspect of this work is how we handle 
black hole singularities. Most existing approaches ei- 
ther use some form of boundary condition on a trapped 
surface on or inside each black hole (see for exam- 
ple [40, 50, 75]), or use a slice that maps the inte- 
rior region of the computational domain for each black 
hole to either part (so called "trumpets" [76]) or all 
("punctures" [77]) of a different asymptotically flat re- 
gion spanned by an Einstein-Rosen Bridge (for a novel 
variant that does not require separation of the metric into 
a background piece and conformal factor see [78]). Here 
we follow an alternative approach where some distance 
inside the apparent horizon of each black hole we replace 
the vacuum interior with an (unphysical) distribution of 
stress-energy to regularize the interior metric. This is 
similar to the "turduckening" evolution scheme [79, 80]. 
However, since we use excision to subsequently evolve 
the initial data, with the excision surface chosen to en- 
tirely contain the unphysical matter, here it is merely a 
device to set up a simple initial data problem without ex- 
plicit interior boundary conditions or singularities. Note 
however that if we were to solve the ID on a domain 
with traditional excision surfaces inside each black hole, 
we would (assuming a well-posed elliptic problem) ob- 
tain the same solution exterior with appropriate excision 
boundary conditions, though the mapping between some 
unphysical interior and appropriate boundary conditions 
would be non-trivial and in general non-unique. 

An outline of the rest of the paper is as follows. In 
Sec. II we review the CTS formulation, describe our 
method for choosing the metric and fluid free data, out- 
line the scheme for regularizing black hole solutions, and 
describe how we numerically solve the constraint equa- 
tions using a multigrid solver. In Sec. Ill we present ex- 
amples of initial data obtained with our solver for quasi- 
circular, eccentric, and ultrarelativistic mergers of com- 
pact objects. Finally, we comment on our results and 
discuss possible future improvements in Sec. IV. In the 
Appendix we give some details on how we treat mesh re- 
finement boundaries in our multigrid algorithm. We use 
geometric units where Newton's constant G = 1 and the 
speed of light c = 1. 



II. COMPUTATIONAL METHODOLOGY 

A. Conformal thin-sandwich equations 

To formulate the initial data problem for general rel- 
ativity, we start by foliating spacetime with a family of 
spacelike hypersurfaces T, t parametrized by t. The nor- 
mal vector to these surfaces nf 1 and the generator of time 
translations t M satisfy 
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where a is the lapse and /3 M is the shift, which is tan- 
gent to S t {n^fi^ = 0). We use the standard convention 



where Greek indices run through {0,1,2,3} and repre- 
sent the full spacetime coordinates, while Latin indices 
run through {1,2,3} and represent coordinates intrinsic 
to a given spatial hypersurface. Using the orthogonal 
projection operator i_ M v = (5 M „ + n. M n„, we obtain the 
induced metric on S t , 7y = g^ v _L M , iy j, where g^ v is 
four-dimensional spacetime metric. The line element can 
be written in terms of these quantities as 
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-a 2 dt 2 + 7^ (dx l + I3 l dt){dx j + ftdt). (2) 



The extrinsic curvature of a slice S t can be written in 
terms of a Lie derivative as 
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Projecting the Einstein equations onto to the hypersur- 
face S 4 one obtains the constraint equations 
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DjK ij - D l K = 8np l 
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where K = j lJ Kij, R and Di are the Ricci scalar and co- 
variant derivative associated with 7y, respectively, and 
E and p % are the energy and momentum density as mea- 
sured by an Eulerian observer, respectively. 

In the language of the 3 + 1 decomposition, initial data 
for the Einstein field equations (and any matter evolu- 
tion equations) are a set of 20 functions representing the 
components of a, j3 l , 7^, -Kjj, E and p l on the initial slice 
S t that together satisfy the constraints (4-5). Though, in 
principle, there are numerous conceivable ways of com- 
ing up with consistent initial data, it is challenging to 
separate freely-specifiable versus constrained degrees of 
freedom in a manner where the underlying physical in- 
terpretation of the free data is transparent, and where the 
choice of the free data leads to a well-posed set of con- 
straint equations. The CTS method [37] is a prescription 
for this separation of degrees of freedom that begins with 
a conformal decomposition of the spatial metric and the 
extrinsic curvature. Introducing the conformal factor \l/ 
we define 
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where j 1 ^ = ^{ffi — 'lkii M ) is defined to be trace- 
less, the overdot indicates a time-derivative, a = ^~ 6 a, 
and R and Di are the Ricci scalar and covariant deriva- 
tive associated with 7^, respectively. With these defini- 
tions we can rewrite (4) and (5) in the CTS form as 

Di&V - + -A ti A lj ^- 7 - — * 5 = -2n^- 3 E 
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DjA ij - -^D l K = 8np l (9) 
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with p i = ^ 10 p\ E = ^ S E. Initial data is obtained by 
solving this system of four elliptic equations for VP and (3 l 
(upon substitution of (7) into (9)), where jij, 7 y , K, a, 
E, and p 1 are the "free data" which can be freely specified 
to reflect the physical system under investigation. 
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B. Superposed free data 



Under the conformal thin-sandwich method one is free 
to choose any values for 7^, 7^, a, K, E, and p l for 
which a solution can be found. In this section we out- 
line our method for determining this free data in order to 
construct initial data representing binary systems. The 
basic idea is as follows. Since solutions to the Einstein 
equations representing isolated compact objects (black 
holes, stars, etc.) are well known, and since if the sepa- 
ration between the objects is not too small the solution 
describing two compact objects is well-approximated by 
superposing the two isolated solutions, we therefore set 
our free data using such a superposed solution and then 
solve the constraint equations in order to obtain the non- 
linear correction. 

There are many ways to combine the metrics repre- 
senting isolated compact objects. The method we use is 
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based on the 3 + 1 splitting. Let j> 
^H 1 ) represent the spatial metric, time derivative of the 
spatial metric, lapse, and shift, respectively, of the first 
isolated solution (e.g. a boosted black hole or neutron 
star solution) and similarly for the second isolated solu- 
tion. Then, we construct the following quantities: 
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where is the flatspace metric. This particular con- 
struction will break down if a( 8up ) < or det[ 7 ^ up ^] < 
anywhere in the domain, which would then require some 
other way of combining the metrics, for example, us- 
ing distance-weighted attenuation functions as in [81]. 
(In [62] it was also found necessary to enforce a desired 
asymptotic fall off of the superposed metric, due to the 
use of a co-rotating frame.) However, these conditions 
are not violated for the cases considered here. From the 
above quantities, we then calculate the free data we will 
use when solving the CTS equations from the usual rela- 
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For initial data with matter we use a similar method. 
We set E and p 1 by superposing the energy and momen- 
tum density of the two objects (we do not consider sit- 
uations where they would both be non-zero at the same 
point). For some cases (in particular for the ultrarela- 
tivistic boosts) , we rescale the momentum density so that 
its magnitude with respect to the superposed metric "/ij 
is equal to the magnitude of the original momentum den- 
sity with respect to the metric of the isolated object (7!- 

or "1^)- This ensures that E 2 and p % pi have the same 
ratio as the isolated objects. This is important since the 
choice of conformal scaling of the energy E = E^> s was 
designed to ensure that if the conformal quantities satisfy 
the dominant energy condition, jijp 1 ^ < E, then so 
will the rescaled quantities following the solution of the 
constraints. 



C. Regularizing black hole solutions 

In cases where black holes are a part of the physical sys- 
tem, the divergences at the black hole's singularity must 
be addressed. As discussed in the introduction, there are 
several ways to deal with this issue in the initial data 
problem; the approach we take here is to explicitly mod- 
ify the metric of an isolated (prior to superposition) black 
hole solution inside the horizon to take a prescribed, reg- 
ular form. The regularized region will not in general 
satisfy the vacuum constraint equations, and to avoid 
a singular conformal factor and shift vector components 
when solving the constraints with such background data, 
we introduce unphysical energy-momentum in the union 
of black hole interiors so that these regions automatically 
satisfy the constraints, albeit with the unphysical interior 
matter source. 

We start with a single, unboosted spinning black hole 
spacetime in horizon penetrating coordinates (for the re- 
sults described here we use the harmonic form of Kerr 
derived in [82], though we have also tried it using Kcrr- 
Schild coordinates without difficulty), so that the only 
divergences in the metric components are well within the 
horizon. We then choose a surface that encloses the sin- 
gular region, yet is strictly inside of the event horizon. 
The interior of this surface we call the regularization re- 
gion. Outside the regularization region we do not modify 
the metric. Inside, there are many conceivable ways to 
alter the metric to eliminate the divergences. The simple 
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approach we take is to promote the black hole mass M 
and spin a constants to functions of space, and smoothly 
decrease them from their bare values at the regularization 
surface to zero at some surface interior to this. 

Specifically, we introduce a regularization function 



freg(x) 



(6a 



1 x> 1 

- 15a? +10) 1 > x > (18) 
> x 



chosen to be twice continuously diffcrentiable so that the 
consequent unphysical energy is well-behaved. We use a 
Cartesian grid 1 and define A E \x, y, z) = \J x 2 + y 2 + z 2 
as the Euclidean radius for a point with coordinates 
x,y,z. We then replace the mass M and the spin pa- 
rameter a with £(a;, y, z)M and £(a?, y, z)a respectively in 
all the metric components, where 



{r i - E \x,y,z)/r { ^\x,y 1 z)) 



eg 



9out 



(19) 

r + (x,y,z) is the Euclidean radius for the point on the 
event horizon at the same angular direction as (x,y,z), 
<7out defines the outer surface of the regularization region, 
and <?i n the inner surface inside of which the metric is 
Minkowski, with 1 > q out > q in > 0. The shape of 
the regularization region, namely a shrunken form of the 
interior of the event horizon, was motivated by the similar 
volume excised during evolution (though that is based on 
the apparent horizon, and the excision surface is a best- 
fit ellipsoid rather than the exact shape of the apparent 
horizon). The particular values of q on t and q- ln are not 
too important (i.e., give essentially the same solutions), 
the only practical requirements being that q out represents 
a surface within the excision surface we will use during 
evolution, and that q- ln not be too close to g ut, otherwise 
excessive resolution is needed to resolve the transition. 

Once we have an everywhere-regular metric, we super- 
pose it with any other COs to construct the free data as 
described in Sec. II B. We then compute the unphysical 
energy and momentum we will add to the regularization 
regions simply by evaluating (4) and (5) with the back- 
ground, superposed data: 
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-E unp hy S and p unphsy are then added to E and p 1 within the 



1 Note that in the harmonic coordinates of [82] only the region 
with rji > M is represented on the Cartesian grid r > 0, where 
rpc is the radial coordinate of the the metric in ingoing null Kerr 
form. Hence the physical singularity is not on the grid, however, 
the metric components are discontinuous &tx = y = z = 
(tk = A^)> hence regularization is still required. 



regularization regions 2 , and we can then solve the CTS 
equations as usual without any additional special treat- 
ment of these regions. During evolution, we choose black 
hole excision surfaces that entirely contain the regulariza- 
tion regions and unphysical matter. Thus, one can think 
of the unphysical matter as serving as a proxy for what 
would otherwise be boundary conditions for ^ and f3 l 
on excision surfaces. Given a solution to the constraints 
with regularized interiors it is trivial to read off what the 
equivalent (Dirichlet) boundary conditions would have 
been, though the inverse problem of mapping some set 
of desired boundary conditions to interior sources is less 
trivial, and likely not well-posed in general. 



D. Fluid Solutions 

For the applications with (physical) matter considered 
here we use Tolman-Oppenheimer-Volkov (TOV) star so- 
lutions in isotropic coordinates to construct the metric 
free data quantities as well as E and p l . Such solutions 
are derived by assuming a relationship between the pres- 
sure and density P(p), e.g. as given by a polytropic con- 
dition. Once the constraint equations have been solved 
and E and p l found, we determine the new density and 
pressure profiles using this same relationship and solving 
the equation 



(E + P(p))(E-p)-p i p i = 



(22) 



for p, which follows directly from the expressions for the 
energy and momentum density of a perfect fluid. For the 
applications considered here, we do not explicitly impose 
any additional constraints on the fluid quantities (e.g. 
that the fluid be in hydrostatic equilibrium). We leave 
that to future extensions. 



E. Multigrid elliptic solver 

In order to numerically solve the CTS equations we 
discretize (8) and (9) using standard second-order finite 
difference operators and solve them using a full approx- 
imation storage implementation of the multigrid algo- 
rithm with Adaptive Mesh Refinement (AMR) as de- 
scribed in [83]. A multigrid algorithm is characterized 
by a smoothing operation, and a choice of restriction and 
prolongation operators. We use Newton-Gauss-Seidel re- 
laxation for smoothing, and half-weight restriction and 
linear interpolation for the restriction and prolongation 
operators, respectively. These latter operators require 
special treatment on mesh refinement boundaries which 
we outline in the Appendix. Unlike in the evolution 



2 It is also possible to calculate the unphysical energy-momentum 
before the superposition, and add the free data E and p after- 
wards; either approach gives similar results. 
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code, we do not use a compactificd coordinate system 
extending to spatial infinity, rather the initial data nu- 
merical grid extends to a large but finite radius. At the 
outer boundaries we impose boundary conditions that 
* = 1 and P l = /3'( su p). Any points outside of this 
domain on the evolution grid are initialized via extrap- 
olation, assuming a leading order 1/r approach to an 
asymptotically-flat spacetime. 

For some applications we wish to solve for initial data 
with axisymmetry. In order to efficiently solve the con- 
straint equations in these situations we have implemented 
a modified version of the Cartoon method [84] similar to 
that used in [85] . Letting the y-axis be the axis of sym- 
metry, we restrict our computational domain to a subset 
of the half plane (x, y) G (— oo, oo) x [0, oo). We use the 
existence of an axisymmetric Killing vector to express 
derivatives in the z direction in terms of derivatives in the 
x and y directions. On the y-axis we impose regularity 
which gives the following conditions for the constrained 
variables: d y ^ = and d y P x = d y P y = (3 Z = 0. 
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III. APPLICATIONS 
A. Quasi-circular binary black holes 

As a first application of our technique, we generate 
and evolve ID for the (approximate) quasi-circular in- 
spiral of two non-spinning, equal-mass black holes. Our 
present method for providing free data is not designed to 
easily give initial data for quasi-circular inspiral (though 
presumably with sufficient fine-tuning of the boost vec- 
tors this could be achieved), and this basic example is 
mainly to provide a relatively low eccentricity binary, a 
couple of orbits before merger, for comparison to past 
studies. Specifically, we are interested in seeing how the 
bare, input parameters of the black holes changing follow- 
ing solution of the constraints, and how much "spurious" 
gravitational radiation is present in the initial data. 

For the initial data, we use free data set by superpos- 
ing two boosted non-spinning equal mass black holes at a 
coordinate separation of 10M, where M is the sum of the 
isolated black hole masses (which in general will be dif- 
ferent from the irreducible masses of the black holes once 
the constraint equations are solved) . The black holes are 
given purely tangential boost velocities chosen so that, 
when evolved, the black holes undergo a few orbits with 
monotonically decreasing proper separation. The initial 
data grid extends to ±2048M in all three directions. For 
convergence studies of the initial data solver, we use three 
base grid sizes of 33 3 , 65 3 and 129 3 and 12 levels of 
mesh refinement with identical grid structures in each 
case. As expected, the conformal factor and shift vec- 
tor exhibit second-order convergence as shown in Figs. 1 
and 2, as does the residual of (4) and (5). For evolu- 
tion, we use the highest resolution initial data. The 
ID is evolved using the generalized harmonic formulation 
of the field equations, choosing harmonic coordinates at 



FIG. 1: The conformal factor * from BH-BH ID. Upper: * 
on the x-axis which lies on the orbital plane and goes through 
the centers of the black holes. Lower: Differences in $ with 
resolution on the rr-axis, scaled assuming second-order con- 
vergence. 




x(M) 



FIG. 2: Differences in the shift component j3 v with resolution 
on the IE-axis from BH-BH ID, scaled assuming second-order 
convergence. 



t = and transitioning to a damped harmonic gauge 
as described in [86]. The eccentricity is estimated to be 
e « 0.05 based on the evolution of the coordinate dis- 
tance between the centers of the apparent horizons as 
shown in Fig. 3. Though the orbital eccentricity could 
presumably be reduced further by tuning the initial ve- 
locities using methods such as the one proposed in [87], 
we did not attempt to do so for this basic comparison. 
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FIG. 3: Coordinate separation of the centers of the two black 
holes fitted to a function (A-B(t~t Q )) 1/4 +C cos(w(t-to)+0). 
This function combines the decaying orbit due to quadrupole 
radiation with the effects of eccentricity, given by e = C/d(t = 
to) ~ 0.05. Due to early-time gauge effects (a transition from 
harmonic to damped harmonic gauge) we exclude the first 
t = 40M irr from the fit. 



Because of corrections from solving the constraints, the 
sum of the masses of the isolated black holes whose space- 
times we superpose M is different from the sum of irre- 
ducible masses computed from their apparent horizons at 
the beginning of the evolution M; rr . For this particular 
case Mi„/M = 1.21. The ID is constructed using free 
data with non-spinning black holes, and the initial spin 
calculated from the apparent horizons is zero to within 
truncation error (|SyM| H | < 6 x 10~ 3 ). The ratio of the 
irreducible mass of the final black hole after the merger to 
sum of the irreducible masses of the initial black holes is 
Mi TI t/Mi TT — 0.885, and the dimensionless spin param- 
eter of the final black hole is af/Mf = 0.678. Both of 
these values are in good agreement (considering the mild 
initial eccentricity here) with the high accuracy results of 
0.88433 and 0.68646, respectively, from [88]. In Fig. 4 we 
show the gravitational waves from the BH-BH merger. 
The initial spurious part of the signal is of comparable 
magnitude to other ID approaches that do not attempt 
to include gravitational waves from the prior inspiral; see 
for example [70]. 



B. Eccentric compact object mergers 

As another application of this technique, we consider 
constructing initial data describing a dynamical capture 
BH-NS binary. We set the free data using a boosted har- 
monic black hole solution and a neutron star with the HB 
equation of state [89]. Let M be the sum of the masses 
of the isolated black hole and neutron star. We construct 
initial data for a 4:1 BH-NS binary by setting the boost 



FIG. 4: The real and imaginary components of the / = 2, 
m = 2 spin- weight —2 spherical-harmonic of r\E , 4 extracted 
at a radius of 105Af. Time is measured from the beginning 
of the simulation. 



velocities to correspond to a Newtonian orbit with eccen- 
tricity e = 1 and periapse distance r p = 5M at various 
initial separations d. We keep the mass and spin that we 
use for the black hole component of the free data fixed 
at 0.8M and — 0.4M, respectively, (where the negative 
sign indicates that the spin is retrograde with respect 
to the orbital angular momentum) and the mass of the 
neutron star component of the free data fixed at 0.2M. 
The spin and masses will receive corrections from solving 
the constraint equations and with decreasing d these will 
differ more and more from the input parameters of the 
free data. The input parameters can of course be tuned 
to achieve desired values in the final solution. However, 
since here we are mainly interested in quantifying this 
difference, we keep them fixed. We use a grid extending 
from -1600M to 1600A/ in each dimension where the base 
level is covered by 257 3 points and there are 9 additional 
levels of mesh refinement, each with a refinement ratio of 
two. We solve for data with initial separations d/M = 
15, 25, and 50. In Table I we show the maximum dif- 
ference of the conformal factor from unity as well as the 
actual ADM (Arnowitt-Deser-Misner) mass, black hole 
mass and spin, neutron star rest-mass, and induced neu- 
tron star density oscillations for these three different sep- 
arations. We can see that even at a separation of 15 M 
the difference between input and final parameters is small 
— at the level of a few percent. At such separations, how- 
ever, the oscillations induced in the neutron star by the 
initial setup become large. This problem could be reme- 
died by adding additional constraints to the matter, for 
example requiring it to satisfy an equilibrium version of 
the Euler equations. 

We evolve the initial data past merger using the same 
methods, gauge, and three different resolutions as in [15]. 
Unless otherwise stated all quantities are from the high 
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FIG. 5: The L 2 -norm of the constraint violation, C a = H a — 
Ox a , in units of 1/M for the d — 15M BH-NS merger in 
the 100M x lOOM-region around the center of mass in the 

equatorial plane (i.e. \J J \\C a \\ 2 d 2 x/ f d?x ). This is shown 
for low, medium, and high resolutions where the latter two 
are scaled assuming second-order convergence. 



FIG. 6: The log of the magnitude of the I = 2, m = 2 spin- 
weight —2 spherical harmonic of r^4 for BH-NS simulations 
with different initial separations d. The value of ^4 was ex- 
tracted on a sphere of radius 100M and is shown starting at 
the beginning of the simulation and continuing past merger. 
The waveforms have been aligned so that the peaks occurs at 
time 0. 



resolution runs. In Fig. 5 we show the norm of the con- 
straints throughout the evolution of the d = 15M ID 
at the different resolutions. The single highest, resolu- 
tion ID is used for all evolution runs, so the fact evolu- 
tion constraints are converging to zero indicates that the 
truncation error of the ID is at least as small as that of 
the highest resolution evolution. In Fig. 6 we plot the 
amplitude of the gravitational waves measured from the 
three different evolutions to show the amount of spuri- 
ous gravitational radiation this method of constructing 
ID introduces. The level of spurious gravitational radia- 
tion decreases with increasing separation and in all three 
cases is small — an order of magnitude or more below 
the physical signal of interest. After the passage of the 
spurious gravitational radiation, the gravitational wave 
signal from all three initial separations is approximately 
the same, though there are small differences due to the 
changes in parameters indicated in Table I, and because 
we are starting the systems along different points of a 
Newtonian trajectory. 



C. Ultrarelativistic initial data 

As a final application, wc consider the problem of 
specifying ID for ultrarelativistic collisions. The study 
of the collision of objects where kinetic energy domi- 
nates the dynamics of the spacetime is of considerable 
interest to super-Planck scale particle collisions, as ar- 
guments suggest classical Einstein gravity will be ade- 
quate to describe the process [90-92]. The hoop con- 



jecture [93] predicts that the generic outcome of a suffi- 
ciently ultrarelativistic collision will be black hole forma- 
tion, and this, together with suggestions of a TeV Planck 
scale [94, 95], imply that, if such a scenario describes 
nature, the Large Hadron Collider (LHC) or cosmic ray 
collisions with Earth could produce black holes [96-98]. 
Though to date no signs of black hole production have 
been observed [99, 100], the nature of the kinetic energy 
dominated regime in general relativity is of interest in its 
own right, and has largely been unexplored. 

Initial data describing such systems will be far from 
equilibrium and one cannot assume that the solution is 
time-symmetric or quasi-static. It is instructive to recall 
the Aichclburg-Sexl [101] solution describing a gravita- 
tional Shockwave. The solution can be obtained from a 
boosted Schwarzschild solution by simultaneously taking 
the mass to zero and the boost parameter to infinity, 
while keeping their product constant and finite. Two 
such oppositely boosted solutions can be superposed to 
obtain a new solution that is valid up until collision. 
Though it is not clear how applicable this is to the non- 
limiting case, this suggests that superposition may be a 
good approximation to describing such spacetimcs. 

Here we consider the specific example of the setup for 
a head-on collision of two fluid star solutions. We use 
the method described in Sec. II B to construct free data 
from two r = 2 polytropic TOV star solutions that have 
unboosted mass M* and a compactness (ratio of mass-to- 
radius) of C = 0.01. The stars are boosted towards each 
other with boost factor 7 = 10. We consider a sequence of 
solutions at various initial coordinate separations d. We 
take advantage of the axisymmetry of the problem and 
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d/M 


max(|* - 1|; 


M /M 0>oo 


A/bh/M a BH /M A/adm/M p OBcm . (%) 


15 


0.0155 


1.077 


0.832 -0.398 1.051 14.3 


25 


0.0092 


1.049 


0.818 -0.402 1.030 9.0 


50 


0.0046 


1.028 


0.808 -0.399 1.017 4.5 



TABLE I: Characteristics of BH-NS initial data with Newtonian orbital parameters r v = 5M and e = 1 with three initial 
coordinate separations d. Here max(|\I/ — lj) is the maximum deviation over the entire domain of the conformal factor from the 
background free-data value of unity, Mo/Mo j00 is the rest-mass of the neutron star compared to its isolated rest-mass, Mbh/M 
and aBH/M are the black hole mass and spin parameters measured from the apparent horizon relative to the initial total mass 
M of the free data, Madm/M is the relative ADM mass of the solution, and Poscm. is the relative magnitude of the oscillation 
in time of the maximum rest-mass density of the neutron star induced by the ID construction. 



use [-2000M, 2000M] x [0,2000M] where M = 
as our computational domain. The base level is covered 
by 1025 x 513 points and there are 9 additional levels of 
mesh refinement. To test convergence we also consider 
two lower resolutions with grid spacing 2 and 4/3 times 
as coarse. 

Using the method for specifying free data described in 
Sec. II B, as d — > oo we expect the corrections from solv- 
ing the constraints will go to zero: ^ — > 1, the magnitude 
of the coordinate velocities of the stars \v\ will approach 
y/l -7~ 2 , the ADM mass Madm will approach M, and 
the total rest-mass Mq will approach the sum of the rest- 
masses of the isolated stars JWo i00 . In Fig. 7 we show 
how all these quantities change with coordinate separa- 
tion. We can see that it is possible to solve for ID where 
the stars are quite close together, though the corrections 
become large, and in particular the ADM mass decreases 
quite significantly. 

To give an indication of the numerical errors on these 
quantities we can compare the values obtained at the 
highest resolution to the Richardson extrapolated val- 
ues using all the three resolutions. For example, for the 
smallest separation d = 1.56M, we have maxd^ — 1|) = 
0.05326 (0.05325) and max|v| = 0.530149 (0.530153) 
where the values in parentheses are the Richardson ex- 
trapolated quantities (which are consistent with approx- 
imately second-order convergence). 

We can compare the above method of constructing free 
data for this case to a conformally flat method. Specif- 
ically, we set all the metric free data quantities to their 
flatspace values, and set E and p l for each star to a 
special-relativistically boosted density and pressure pro- 
file taken from the TOV solution. In Fig. 8 we show the 
same quantities as in Fig. 7 but using this conformally 
flat method. In this case the corrections from solving 
the constraints will not go to zero with infinite separa- 
tion since all the non-trivial geometry is coming from the 
conformal factor. Hence the energy-momentum will be 
substantially rescaled at any separation. Also in contrast 
to the first method, the maximum of \*f? — 1| occurs for 
\& > 1 instead of $ < 1, which means E and p l will be 
smaller than their conformal counterparts. With confor- 
mally flat ID it is also possible to solve for stars close 
together though, as in the previous method, the ADM 



mass decreases steeply. It should also be noted that be- 
cause of the large shift vector obtained with the second 
method, the coordinate velocity is substantially greater 
than one, which may make it more challenging to numer- 
ically evolve. 

A full characterization of this ultrarelativistic collision 
ID will require evolution, which we will present in up- 
coming work [102]. 



IV. CONCLUSIONS 

We have outlined a general method for constructing 
initial data based on superposition and the CTS for- 
mulation of the constraint equations, and demonstrated 
the method with some example solutions. Though there 
are numerous existing applications of the CTS method, 
and superposition has been proposed before, some of the 
novel aspects of the work presented here include adding 
neutron stars to the prescription, regularizing the interi- 
ors of black holes with (unphysical) matter sources, and 
applying it to regimes not yet studied before, namely 
initial data for high-eccentricity binary mergers and ul- 
trarelativistic collisions. For astrophysically relevant bi- 
naries we find that superposition of single, isolated com- 
pact object solutions works well in the sense that non- 
linear correction from solving the constraints are rela- 
tively small for larger initial separations, implying that 
superposition is a good start to attain more astrophysi- 
cally realistic initial data (for example, by adding prior 
gravitational wave information as in [70] to the super- 
posed background data for quasi-circular or low eccen- 
tricity inspirals). Including neutron stars, we find that 
the superposition effectively induces oscillations in the 
stars. This again is small for large separations and hence 
a good approximation to dynamical capture binaries. 
However, practical application to low eccentricity inspi- 
rals will likely require that the CTS equations be supple- 
mented with some form of quasi-equilibrium equations 
for the hydrodynamics (as in many existing ID methods, 
for example [54, 59, 72]). 

For the ultrarelativistic boost examples we are able to 
obtain solutions to the CTS equations with superposed 
and conformally flat data well into the kinetic energy 
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to the conformal constraint equations beyond constant 
mean curvature slicing [103] (and in some cases, such as 
the extended CTS equations [104] there are known exam- 
ples of non-uniqueness [105]), it is interesting that we are 
able to obtain solutions in this highly non-linear regime. 
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FIG. 7: Various quantities from ultrarelativistic collision ID 
with 7 = 10 made using the superposition method for con- 
structing free data. From top to bottom the quantities shown 
are the maximum (over the entire domain) difference of the 
conformal factor from unity, the maximum coordinate veloc- 
ity of the fluid, the total rest-mass, and the ADM mass. All 
quantities are shown as a function of d, the coordinate separa- 
tion between the two stars. For all these cases the maximum 
of |$ — 1| occurs for values of >& that are less than unity. One 
might expect these quantities to approach their infinite sep- 
aration limits as 1/d for large d; the dotted lines show such 
1/d curves for each quantity matched to the d/M = 50 point. 
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dominated regime (7 = 10) for sufficiently large initial 
separations. At smaller separations we are still able to 
obtain solutions. However, for these initial data sets the 
corrections to the metric and fluid properties become 
large, and it is less clear how to separate the total energy 
of the spacetime into kinetic energy, rest-mass energy, 
etc. This will require evolution to resolve, and we leave 
that to future work. Nevertheless, given that there are 
few results on the uniqueness and existence of solutions 
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FIG. 8: Same as Fig. 7 but with conformally flat data. For 
all these cases the maximum of |$ — 1| occurs for values of 
that are greater than unity. 
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Appendix: Multigrid AMR interpolation 

A multigrid algorithm requires a restriction operator 
to inject quantities from finer to coarser grids as well as 
a prolongation operator to interpolate corrections from 
coarser grids to finer grids (see e.g. [106]). For our multi- 
grid algorithm we use half-weight restriction as our re- 
striction operator. In three dimensions half-weight re- 
striction can be written as 

/hw = fi, j)k + ^ (Afxx + A/„„ + Af zz ) (A.l) 
where 

A xx f — fi+i,j,k - 2/jj. fe + fi-i,j,k (A. 2) 

and similarly for the y and z directions. Note that A xx f 
divided by h 2 (where h is the grid spacing) is a second- 
order approximation for d 2 f. On AMR boundaries where 



the full stencil is not available we must modify the above 
expression. For example on a negative x-boundary we 
replace A xx f by a right-handed second derivative stencil 

A xx f = 2fi,j,k — 5fi-i,j,k + 4/i-2j,fc — fi-3,j,k (A. 3) 

and so on for the other directions. This ensures not only 
that /hw is a second-order representation of /, but also 
that /hw is smooth to 0(/i 4 ) on AMR boundaries. Hence 
if second derivatives of /hw are computed including re- 
stricted boundary points in the stencil the error will be 
0(h 2 ). 

We use linear interpolation as our prolongation oper- 
ator. However, after applying a correction from a coarse 
grid, we reset the values on the AMR boundaries of the 
fine grid for the points that do not exist on the coarse 
level with fourth-order interpolation using those points 
that do. We found this higher order interpolation to be 
beneficial as we do not relax the points on the boundary. 
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